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The Berezinskii-Kosterlitz-Thouless (BKT) transition in two-dimensional planar rotator and XY 
models on a square lattice, diluted by randomly placed vacancies, is studied here using hybrid 
Monte Carlo simulations that combine single spin flip, cluster and over-relaxation techniques. The 
transition temperature T c is determined as a function of vacancy density p vac by calculations of the 
helicity modulus and the by finite-size scaling of the in-plane magnetic susceptibility. The results for 
T c are consistent with those from the much less precise fourth-order cumulant of Binder. T c is found 
to decrease monotonically with increasing p va c, and falls to zero close to the square lattice percolation 
limit, pvac ~ 0.41 . The result is physically reasonable: the long-range orientational order of the 
low-temperature phase cannot be maintained in the absence of sufficient spin interactions across the 
lattice. 

PACS numbers: 75.10.Hk, 75.30.Ds, 75.40.Gb, 75.40.Mg 



INTRODUCTION: SPIN-DILUTED PLANAR 
SPIN MODELS 



It is well known that vortices are fundamental in- 
gredients in the_Berezinskii-Kosterlitz-Thouless (BKT) 
phase transition.mm The simplest models exhibiting this 
transition are the pure planar rotator model (PRM) 
and XY-model. For these models the transition takes, 
place at critical temperatiujes ksTxr/JS 2 = 0.898a 
and k B T KT /JS 2 = 0.699jala respectively (J is the ex- 
change constant, S the spin length). Recently, the 
study of topological excitations such as vortices and soli- 
tons in two-dimensional magnetic lattise^ coi | itaining de- . 
fects has received a lot of attention .IffiuHOE3ll3llil^ 
Such interactions must have interesting consequences for 
the static and dynamical properties of easy-plane mag- 
nets. Analytical and numerical calculations have shown 
that vortices, axe.,aljtracted and pinned by nonmagnetic 
impurities. EjOEjEj In fact, the vortex energy is lowered 
when pinned at a vacancy, resulting in— greater prefer- 
ence of single vortexcSl and vortex-pairtH formation on 
vacancies. Of course, this leads to an overall increase 
of the system disorder. All of these factors conspire to 
reduce the BKT transition temperature with increasing 
vacancy densit y, a s has already been seen in calculations 
from Refs. |lq , [19|j2C| for planar spin models on a two- 
dimensional square lattice (see also analytical results us- 
ing the self-consistent harmonic approximation of Ref. 
^T[ ). The important question here is, at what vacancy 
density is the transition temperature reduced to zero, 
so that the system is always in the high-T disordered 
phase? This would mean a situation in which there is 
no low-temperature phase of long-range orientational or- 
der, characterized by spin-spin correlations decaying as 
a power law with distance, and a finite absolute magne- 
tization (| Si\) in the thermodynamic limit. 



Calculations of the helicityjnodulus for the planar ro- 
tator model by Leonel et alx3 indicated that a critical 
vacancy density p vac = p c rj 0.3 causes the critical tem- 
perature T c to drop to zero. It means that the critical 
temperature would vanish at p c = 1 — p c w 0.7, which is 
above the site percolation threshold, p p t = I — p p t t^-P-59, 
for a planar square lattice. Lozovick and Pomirchi,EHI also 
using the jump in the helicity modulus, have found that 
the BKT phase transition occurs above the percolation 
threshold in a dilute system of Josephson junctions (us* 
ing bond dilution). On the other hand, Berche et alx3 
calculated the decay of the spin-spin correlation function 
and its related exponent, 77, and considered the transi- 
tion temperature to be located by T](T C ) = 1/4. Those 
results suggested that the critical density is closer to 0.41 
(the number associated with the percolation limit for the 
square lattice). The Monte Carlo calculations for this 
problem naturally are particularly difficult, especially be- 
cause the interesting region occurs at very low tempera- 
ture. Furthermore, the statistical fluctuations due to the 
random choice of positions for the vacancies further in- 
creases the numerical noise in the calculations - this effect 
itself becomes particularly troublesome especially when 
Pvac surpasses 0.3 (30%). As such, it seems important 
to make more reliable estimates for the critical vacancy 
density based on improved MC calculations here. 

The calculations mentioned above concern the planar 
rotator model (two-component spins lying in xy plane). 
In a specialized model with repulsive vacancies, Wysin 
calculated the reduction of T c in an easy-plane Heisen- 
berg model, with three-component spins with anisotropic 
couplings of their components.E3 The vacancies were not 
allowed to be on nearest or next nearest neighbor lattice 
sites, which made it possible to calculate the vorticity 
density in the model. However, that calculation did not 
concern itself with the determination of the critical va- 



2 




k R T/JS 

FIG. 1: Application of the fourth order cumulant (^J) for 
estimating T c , for the PRM at 4% vacancy concentration. The 
data we re ob tained using the Monte Carlo approach described 
in Sec. II B. The inset expands the view near the critical 
temperature (k B T c /JS 2 « 0.815). 



cancy density, because the constraint of repulsive vacan- 
cies limits the possible vacancy density to be less than 
18%, well below the critical value. Therefore, for com- 
parison with the planar rotator, we also consider here 
the vacancy effects in the (three-component) XY-model, 
with randomly placed non-repulsive vacancies. 

After describing the model Hamiltonians, we give an 
overview of the different methods used to estimate the 
transition temperature. This is followed by some spe- 
cific comments on the Monte Carlo schemes applied to 
this problem. The data obtained for the planar rotator 
and XY models are presented, followed ultimately by our 
conclusions. 



II. MODEL HAMILTONIANS 

The spin models under consideration can be described 
by the Hamiltonian 



H 



-J ^VjfT,' 

(i,3) 



(s?s x 



s f s j) 



(1) 



where (i, j) indicates nearest neighbor sites of an L x L 
square lattice, and J is the ferromagnetic exchange cou- 
pling between spins Si and Sj. The spins S% have two 
components for the planar rotator model and three com- 
ponents for the XY-model; in the latter case, however, 
only the xy components are coupled. The occupation 
variables a take the values 1 or depending on whether 
the associated site is occupied by a spin or vacant. A frac- 
tion p vac of the sites are chosen randomly to be vacant. 
It is important to realize, however, that the Monte Carlo 
calculations here must make adequate averages over dif- 



ferent choices of the vacancy positions, for a chosen den- 
sity. 

The planar rotator model has effectively a single degree 
of freedom per site - the angle of the spin within the xy 
plane. The main distinction of the XY model is the pres- 
ence of the extra S z components, which act as degrees of 
freedom, but do not appear in the Hamiltonian. The XY 
model therefore involves two degrees of freedom per spin. 
This increases the entropy effects at a given temperature 
and results in a lower T c compared to the planar rotator. 
The MC algorithm for the XY model must involve the 
possibility to change all three spin components for the 
XY model, while preserving the spin length. 



A. Physical properties leading to T c 

The lack of significant sharp peaks in the thermody- 
namic quantities versus temperature T for these models, 
especially in finite Lx L lattice systems, means that pre- 
cisely locating T c is difficult. Therefore, it is useful to 
apply several different approaches, all essentially based 
on the scaling of the thermodynamics with the system 
size or edge length L. 

As the Monte Carlo algorithm proceeds, the to- 
tal system instantaneous in-plane magnetization M — 
(M x , My) is observed, 



M = J2aiSi 



(2) 



Additionally, statistical fluctuations give the susceptibil- 
ity components for temperature T, 



X aa = {{Ml) - (M a ) 2 )/ (NT). 



(3) 



The number of spins in the system is N = (1 — p va , c )L 2 . 
The average of \ xx an d X yy defines the in-plane suscep- 
tibility, 



1 



x = ^(x xx + x vv )- 



(4) 



A rough estimate of T c can be obtained from-the, size- 
dependence of Binder's fourth order cumulantSEil JJl, 
defined by 



U L = 1- 



2<M2 + M2)2 • 



(5) 



For any L, the asymptotic values are Ul(T <C T c ) = 0.5, 
Ul(T 3> T c ) = 0. At the critical temperature, Ul is 
approximately independent of the system size, hence, T c 
can be estimated from the crossing point of curves of 
Ul(T) for various L. An example of such crossing be- 
havior is given in Fig. H for the PRM at p vac = 0.04. 
In practical application, due to the statistical uncertain- 
ties, there is usually no clear crossing point, especially at 
higher vacancy concentrations. Instead, T c is very close 
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FIG. 2: Typical application of the helicity modulus for esti- 
mating T c , for the PRM at 4% vacancy concentration. The 
dashed line is Eq. Q . The inset shows how the crossing points 
occur slightly above the critical temperature (fcsTc/JS' 2 ps 
0.815). Error bars are smaller than the symbols. 



to the point where different curves of Ul (T) begin to sep- 
arate from the low-T asymptotic value. Although very 
reliable, this approach is not very accurate, and requires 
MC calculations for many temperatures near T c . 

Another approach to determine T c is based on the cal- 
culation of the helicity modulus per spin, Y(T). It is a 
measure of the resistance to an infinitesimal spin twist A 
across the system along one coordinate, defined in terms 
of the dimensionless free energy, / = F/(JS 2 ), 



T = 



i a 2 / 



(0) 



NdA 2 

Any general model Hamiltonian leads to the expression, 
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where f3 = (fc^T) -1 is the inverse temperature. For ei- 
ther the planar rotator or XY model, the required oper- 
ators to be averaged (in limit A — > 0) can be expressed 
using the Cartesian spin components, 

G ^0A=Y, Wi ' 4 ) ~ S i S D > ( 8a ) 



(id) 



G c 
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where e-ij is a unit vector pointing from site j to site 
i. The sum determining G s only includes pairs of lat- 
tice sites displaced by ±x. Furthermore, one expects the 
mean of G s to be quite small, while its fluctuations do 



FIG. 3: Log-log plot of susceptibility versus system edge 
length L, for the PRM at 4% vacancy concentration. The 
curves correspond to different values of the dimensionless tem- 
perature r = UbT/JS 2 . Lines are guides to the eye; errors 
are smaller than the symbols. Least squares fits were used to 
determine the slopes, (2 — 1]), producing rj(T) as seen in Fig. 
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contribute to the helicity formula (Q). The sum for G c is 
seen to be proportional to the original Hamiltonian. 

According to renormalization-group theory, the helic- 
ity modulus in an infinite system jumps from the finite 
value (2/7r)fcsT c to zero at the critical temperature. As- 
suming this applies also to the spin-diluted model, as 
argued in Ref. |l8|, T c can be estimated from the intersec- 
tion of T(T) and the straight line, 



-knT. 



(9) 



The trend in the intersection point with increasing L can 
be observed, as shown for the PRM at p vac — 0.04 in 
Fig. Generally speaking, the MC data for T(T) show 
a steeper drop in the critical region as L increases. The 
larger system size used, the lower will be the intersection 
point and estimated T c . Hence this method will always 
lead to an over-estimate of T c . 

Finally a third approach was also applied for estimat- 
ing T c , employing the finite size scaling of the in-plane 
susceptibility, as used in a pure XXZ model by Cuccoli 
et alB and— in the same model with repulsive vacancies 
by WysimEl In the absence of vacancies, it is the most 
precise method, because the statistical fluctuations in \ 
can be reduced by extended MC averaging much more 
effectively than those of the helicity modulus or Binder's 
cumulant. We assume that near and below T c , a power 
law scaling of the susceptibility holds, even in the pres- 
ence of vacancies, 



L 2 -\ 



(10) 



where r\ is the exponent for the in-plane spin correlations 
below T c (see Ref. |). By using this equation with cal- 
culations at several lattice sizes, the exponent r\ can be 
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FIG. 4: Application of the correlation exponent r\ for esti- 
mating T c , for the PRM at 4% vacancy concentration, de- 
rived from using systems of sizes L — 16, 32, 64, 96. The 
inset shows how the critical temperature was estimated as 
k B T c /JS 2 w 0.815. 



fitted as a function of temperature. An indication of how 
X scales with system size is given in Fig. pi again for the 
PRM at 4% vacancy concentration. One can note clearly 
how the exponent (2 — 77) (slope of log-log plot for x(-^)) 
decreases as the temperature increases, especially rapidly 
as T passes the transition temperature. 

For the pure PR and XY models (no vacancies), the 
transition is located at the temperature where rj(T) = 
1/4. Then, under that assumption that the vacancies do 
not change the basic symmetries in the transition, but 
only increase the effective entropy present, we can expect 
that the transition can be located in the same way under 
the presence of vacancies, solving 



1 

4' 



(11) 



In the absence of any particular theory for the model with 
vacancies, this can be expected to be a reasonable defi- 
nition for T c . Analysis of powep4aw fits of the spin-spin 
correlations in the diluted PRMO also demonstrated that 



T c occurs very close to the temperature from Eq. (11). 
Its validity is further verified here by the comparison with 
the results for T c due to the hclicity modulus, and due to 
Binder's cumulant, the latter of which is reliable for any 
kind of model, with or without vacancies. Fig. || shows 
its application for the PRM at 4 % vacancy concentra- 
tion, leading to kgT c /JS 2 w 0.815, exactly consistent 
with the estimate from Binder's cumulant (Fig. [I]). 

We also note, that for the pure PRM (no vacancies), 
this fitting of 77, using systems as large as L = 160, leads 
to the estimate k B T c /JS 2 = 0.907±0.004, slightly higher 
than that from Ref. ^. 



FIG. 5: The helicity modulus for the PRM at 33% vacancy 
concentration for system sizes indicated. The dashed line is 
Eq. §). 
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FIG. 6: Binder's fourth cumulant for the PRM at 33% va- 
cancy concentration for system sizes indicated. UbTc/ JS 2 w 
0.14 as estimated from the point where the data for different 
system sizes separate. The lines are guides to the eye. 



B. Monte Carlo Scheme 

Thermal averages for a given system size and temper- 
ature were obtained using a hybrid MC approach, in- 
cluding Metropolis single-spin moves and over-relaxation 
movesu that can modify all spin comprowpjts, in com- 
bination with Wolff single-cluster movesc§c3 that mod- 
ify only the xy components. These are based on stan- 
dard approaches .for spin models, as developed in many 
references OoSEH Further details, as applied to a sys- 

Using a 



20 



tern with vacancies, can be found in Ref. 
hybrid method including cluster and over-relaxation is 
very important especially at low temperatures, as it very 
effectively reduces the problems associated with critical 
slowing down. 
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The programming used for the XY model also applies 
to the planar rotator model; it is only necessary to set 
the out-of-plane components S z = and then never allow 
them to change for the PRM. Thus it is straightforward 
to study the two models with essentially the same MC 
approach. 

The calculations were performed for a range of system 
sizes, including L = 16,32,64,96, and 160. For a given 
L x L lattice, the number of vacancies placed at random 
locations was N vac — p vac L 2 (spins removed from system 
or equivalently, set to zero length). For larger systems 
or very low vacancy density, the results are nearly in- 
dependent of the particular random choice of vacancy 
positions. In the general case, however, it is necessary 
to average over equivalent systems (same L, p vac ) with 
different particular choices of the vacancy locations. The 
statistical variations in the averages are most significant 
especially as the vacancy density approaches the critical 
value that forces T c to zero. These statistical variations 
also are strongest in the smaller systems; conversely, the 
larger systems tend to average out these fluctuations, all 
the better as their area increases. Therefore we averaged 
over a number N sys copies of the system, with this num- 



ber taken largest at small system size. 



we used N& 



64,32,8,4, for L 



spectively. For larger density, p v 



> 



For p vac < 0.35, 
16,32,64,96, re- 
0.35, we doubled 



these values for N, 
N sys — 4 for L 



syB , and additionally included runs with 
160. 

For thermal equilibration before calculating averages, 
5000 MC steps (MCS)l3 were applied for small systems 
(L < 50) and 10, 000 MCS for large systems. For each 
of the N sys individual realizations of a given L and p vac , 
averages at one temperature were calculated using be- 
tween 20,000 and 80,000 MCS (JV d ata), with the great- 
est number applied to the larger systems. For example, 
calculation for one temperature of a 16 x 16 lattice at 
Pvac — 0.1 involved an average over 64 x 25,000 = 1.28 
million MCS. On the other hand, one temperature of a 
96 x 96 lattice at p vac — 0.36 involved an average over 
8 x 80,000 = 640,000 MCS. Near 0% vacancy density, 
these MC parameters produce insignificant error bars; 
when p vac exceeds 30%, on the other hand, the resulting 
error bars are considerably greater and resist reduction. 
As suggested above, the error bars in T, \, and Ul can 
be reduced more readily by increasing N sys than by in- 
creasing iVdata when significant vacancy density is present 
(especially at p vac > 0.3). 
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FIG. 7: The helicity modulus for the PRM at vacancy concen- 
trations p V ac indicated in the legend. The dashed line is Eq. 
([]). Part (a) shows the overall trend; error bars are smaller 
than the symbols. Part (b) displays the behavior as the tran- 
sition is extinguished at the critical vacancy concentration. 



nificant. Even so, looking at the trends in the data with 
system size, in the following we show the MC evidence 
that the transition temperature is reduced to zero when 
the vacancy concentration is approximately 40%, for both 
models. 



Planar rotator model 



III. MONTE CARLO DATA 

Calculations were carried for a range of vacancy den- 
sities from zero to 50%. We especially concentrated on 
the region 0.30 < p vac < 0.40, which required the most 
careful analysis. For vacancy density less than 30%, it 
is clear that there is a transition at a finite temperature, 
for both the PR and XY models. At the higher vacancy 
concentrations, statistical errors were generally more sig- 



At low vacancy concentrations (p vac < 0.20), MC re- 
sults for Ul, T, x, an d v(T) bear a great resemblance to 
those shown above for 4% vacancies, with fairly smooth 
dependencies on temperature. The primary modification 
is the general trend of important features towards lower 
temperature with increasing p vac . At higher concentra- 
tions, errors become more significant. For example, the 
helicity modulus at 33% vacancies and various system 
sizes is shown in Fig. |^. In addition to larger relative er- 
rors, the absolute magnitude of T is drastically reduced. 
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FIG. 8: Application of the correlation exponent r\ for estimat- 
ing T c , for the PRM at the vacancy concentrations indicated 
in the legend, derived from scaling of \ f° r systems of sizes 
L — 16, 32, 64, 96. Part (a) gives the rough overall trend and 
part (b) shows how r\ does not fall to the value 1/4 at vacancy 
concentrations greater than 41%. 



It is very clear, however, that the BKT transition is still 
present at this concentration, with ksTc/JS 2 w 0.20 as 
estimated from the crossing point of the L = 96 data. 
This is additionally supported by the corresponding be- 
havior of Binder's cumulant, seen in Fig. [| which gives 
the estimate ksTc/JS 2 « 0.14, somewhat lower, as can 
be expected. 

An indication of the tendency for reduction of T c with 
vacancy concentration is given in Fig. j?], showing a collec- 
tion of results all for L = 96 systems. While these cross- 
ing points consistently overestimate T c , a better view 
of this critical point reduction is provided by the var- 
ious graphs of r/(T) at different concentrations, Fig. [|. 
One can see clearly that once the vacancy concentration 
passes a value around 41%, the fitted value of r\ does not 
fall below the value 1/4, at least for the lowest tempera- 
tures used (k B T/JS 2 = 0.01). 

In Fig. ^, the critical temperatures extracted from 
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(0) and from the crossing of T(T) with Eq. (|) for L = 96 
systems. The inset shows T c approaches the critical 

region. 



TABLE I: Dependence of dimensionless critical temperature 
t c = ksTc/JS 2 on pvac, as derived from the scaling of in-plane 
susceptibility x> an d using rj(T c ) = 1/4. 



pvac 


r (PR model) 


r (XY model) 


0.0 


0.907 ± 0.004 


0.700 ± 0.005 


0.04 


0.815 ±0.005 


0.637 ± 0.005 


0.10 


0.683 ± 0.004 


0.547 ± 0.005 


0.16 




0.453 ± 0.005 


0.20 


0.456 ± 0.004 


0.384 ± 0.005 


0.30 


0.230 ± 0.004 


0.208 ± 0.005 


0.33 


0.153 ±0.007 


0.147 ±0.005 


0.36 


0.093 ± 0.005 


0.087 ± 0.005 


0.38 


0.050 ± 0.006 


0.049 ± 0.005 


0.39 


0.034 ± 0.006 


0.041 ± 0.005 


0.40 


0.019 ± 0.009 


0.018 ±0.007 


0.41 


0.005 ± 0.009 


0.003 ± 0.007 


0.42 


0.0 ± 0.005 


0.0 ±0.005 


0.44 


0.0 ± 0.005 


0.0 ± 0.005 



r\ and Eq. (|Tl]) and from the helicity modulus (using 
L = 96) are shown as functions of vacancy concentra- 
tion. As mentioned above, the crossing point of T(T) 
with Eq. (|) for any finite system always overestimates 
T c , with less error as larger systems are used. The fitting 
of 77 is more reliable; the data shown here used the scal- 
ing fitting of x using systems with sizes L = 16, 32, 64, 96. 
The numerical values of T c estimated using i](T c ) = 1/4 
are summarized in Table Q. The results give strong evi- 
dence for extinction of the BKT transition at a vacancy 
concentration close to 41%. 



7 




0.1 0.2 0.3 0.4 0.5 

k n T/JS 2 

FIG. 10: The helicity modulus for the XY model at 40% 
vacancy concentration for system sizes indicated. The dashed 
line is Eq. (§). 



B. XY model 

The general trends in MC data for the XY model are 
rather similar to those found for the planar rotator. The 
most obvious distinction, however, is that the extra en- 
tropy due to the out-of-plane spin component forces the 
transition temperature to be lower in the XY model, no 
matter what vacancy concentration is considered. 

It is interesting to show some data at 40% vacancy 
concentration, where the transition is seen to occur very 
slightly above zero temperature. In Fig. [l(] the helicity 
modulus for system sizes from L = 16 to L = 160 is dis- 
played. As the data for increasing system size is seen to 
systematically fall to lower values, this graph alone can- 
not undeniably prove the presence of a transition. How- 
ever, when taken in conjunction with the fits for 77, which 
passes the value 1/4 around ksT/JS 2 ~ 0.018, we can 
say that even at 40% vacancy density there occurs a tran- 
sition at finite temperature. This can be seen in Fig. [Il], 
where rj{T) is shown for the various vacancy concentra- 
tions studied. On the other hand, performing the MC 
calculations at temperatures as low as ksT/JS 2 = 0.01, 
the exponent 77 does not acquire such a low value as 1/4 
even for 41% vacancy concentration. 

In Fig. pi t he critical temperatures extracted from 
77 and Eq. (|ll|) and from the helicity modulus (using 
L = 96) are shown as functions of vacancy concentra- 
tion. The numerical values as derived using rj{T c ) = 1/4 
are given in Table ffl. Just as in the PR model, these re- 
sults demonstrate the extinction of the BKT transition 
at a vacancy concentration close to 41%. As the transi- 
tion is controlled by the in-plane spin components, the 
presence of the extra S z component in the XY model 
changes the overall scale of transition temperatures, but 
does not affect the critical vacancy concentration. 

Finally, we can also make some comparison to the XY 
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FIG. 11: Application of the correlation exponent r\ for esti- 
mating T c , for the XY model at the vacancy concentrations 
indicated in the legend, derived from scaling of \ f° r systems 
of sizes L — 16, 32, 64, 96. Part (a) gives the rough overall 
trend and part (b) shows how 77 does not fall to the value 1/4 
at vacancy concentrations greater than 41%. 



model using repulsive vacancies studied in Ref. |2(J. Con- 
siderable data was presented there for the case of 16% 
vacancies. Therefore it is interesting to note how the 
transition temperature is changed if the vacancies are al- 
lowed to be at completely random positions in the current 
model. 

A graph of rj(T) for this case is given in Fig. O, show- 
ing clearly the transition occurring at fc_eT c / JS 2 ^ 0.453. 
Alternatively, and even with less computational effort, 
the transition can be found as done in Refs. p|pc| by plot- 
ting x/L( 2 ~ v ^ vs. T, taking 77 = 1/4, and looking for the 
common crossing point of data at various system sizes. 
This is seen in Fig. |l4|, which gives the same estimate 
for T c . In the repulsive vacancy model at the same va- 
cancy concentration, the transition occurs at a slightly 
higher temperature, k B T c /JS 2 « 0.478 ± 0.001. The re- 
sult is reasonable; there is greater disorder in the model 
with fully random vacancies, hence, requiring less ther- 
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FIG. 12: The critical temperatures versus vacancy concen- 
tration for the XY model, extracted from fits of r\ together 
with Eq. ( p"l] ) and from the crossing of T(T) with Eq. (g) for 
L = 96 systems. The inset shows T c £LS pvac approches the 
critical region. 



FIG. 14: Application of the finite-size scaling of in-plane sus- 
ceptibility x to estimate ksTc/JS 2 w 0.453 (common cross- 
ing point of the data) at 16% vacancy density in the XY 
model, with fully randomly placed vacancies, using exponent 

TJ=l/4. 
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FIG. 13: Application of the correlation exponent r\ for esti- 
mating T c , for the XY model at 16% vacancy concentration, 
derived from using systems of sizes L — 16, 32, 64, 96. The 
inset shows how the critical temperature was estimated as 
k B T c /JS 2 « 0.453. 



mal disordering due to temperature to reach the high- 
temperature phase, (alternatively, the repulsive vacancy 
model has more built-in order and hence requires greater 
thermal energy per spin to reach the high-temperature 
phase.) 



IV. CONCLUSIONS 

Hybrid MC calculations applied to the planar rota- 
tor and XY models on a 2D square lattice show that 
the BKT transition is extinguished (T c — > 0) at a va- 



cancy concentration close to 41%, a number related to 
the percolation limit. Then, although the BKT phase 
transition has an unusual nature, in which the topolog- 
ical long-range order is destroyed by the unbinding of 
vortices, the percolation problem of systems exhibiting 
such a transition must have some similarities to the tra- 
ditional 2D Ising model. In general, the transition tem- 
peratures for the XY model are lower than those for the 
PR model, due to the extra entropy of out-of-plane spin 
motions, but otherwise, the static properties are closely 
related. The transition temperatures were determined 
most precisely using the finite-size scaling of the in-plane 
magnetic susceptibility, under the assumption that the 
spin-correlation exponent rj goes to the universal value 
1/4 at the transition, regardless of the vacancy concen- 
tration. This is equivalent to saying that the presence 
of spin vacancies does not change any fundamental sym- 
metries of the problem. T c calculated this way is com- 
pletely consistent with the corresponding results from 
the helicity modulus and Binder's fourth order cumulant. 
At vacancy concentration higher than 41%, the intrin- 
sic disorder of the system always produces a phase with 
short range correlations that decay exponentially, i.e., the 
usual "high-temperature" BKT phase whose properties 
are strongly determined by the presence of unbound vor- 
tices and antivortices. The lack of percolation across the 
system at p vac > 0.41 disrupts the ability to generate 
topological long-range correlations. It then becomes im- 
possible to lower the temperature adequately to reach the 
ordered phase of very low vortex density, dominated by 
spin waves. 



9 



Acknowledgments versidade Federal de Vicpsa, Brazil. ARP thanks CNPq 

for financial support. 

GMW is very grateful for support from FAPEMIG as 
a visiting researcher and for the hospitality of the Uni- 



Electronic address: wysin@phys.ksu.edu ; URL: tittp:// 
wiw.phys.ksu.edu/~wysiri; Permanent address: Depart- 
ment of Physics, Kansas State University, Manhattan, KS 
66506-2601 

Electronic address: apcrcira® ufv . br 



Electronic address: 3idiney@fisica.ufjf.b1 

V.L. Berezinskii, Sov. Phys. JEPT 32, 493 (1970). 

V.L. Berezinskii, Sov. Phys. JEPT 34, 610 (1972). 

J.M. Kosterlitz and D.J. Thouless, J. Phys. C 6, 1181 

(1973). 

J. Tobochnik and G.V. Chester, Phys. Rev. B 20, 3761 
(1979). 

A. Cuccoli, V. Tognetti, and R. Vaia, Phys. Rev. B 52, 
10221 (1995). 

H.G. Evertz and DP. Landau, Phys. Rev. B 54,12302 
(1996). 

K. Subbaraman, C.E. Zaspel, J.E. Drumheller, Phys. Rev. 
Lett. 80, 2201 (1998). 

C.E. Zaspel, CM. McKennan, and SR. Snaric, Phys. 
Rev.B 53, 11317 (1996). 

M.M. Bogdan and C.E. Zaspel, Phys. Status Solidi A 189, 
983 (2002). 

AR. Pereira, A.S.T. Pires, J. Magn. Magn. Mater. 257, 
290 (2003). 

AR. Pereira, L.A.S. M61, S.A. Leonel, P.Z. Coura, and 

B. V. Costa, Phys. Rev. B 68, 132409 (2003). 
CM. Wysin, Phys. Rev. B 68, 184411 (2003). 

F.M. Paula, AR. Pereira, L.A.S. M61, Phys. Lett.A 329, 
155 (2004). 

AR. Pereira, S.A. Leonel, P.Z. Coura, and B.V. Costa, 
Phys. Rev.B 71, 014403 (2005). 

H. Karatsuji and H. Yabu, J. Phys. A 29, 6505 (1996). 



16 L.A.S. M61, AR. Pereira, W.A. Moura-Melo, Phys. Rev. 
B 67, 132403 (2003). 

17 AR. Pereira, J. Mag. Magn. Mater. 279,396 (2004). 

18 S.A. Leonel, P.Z. Coura, AR. Pereira, L.A.S. M61, and 
B.V. Costa, Phys. Rev. B67104426 (2003). 

19 B. Berche, A.I. Farinas- Sanches, Y. Holovatch, and R. 
Paredes, Eur. Phys. J. B. 36,91 (2003) 

20 CM. Wysin, to appear in Phys. Rev. B (2005). 

21 L.M. Castro, A.S.T. Pires, and J.A. Plascak, J. Mag. Mag. 
Mater. 248,62 (2002). 

22 Y.E. Lozovick and L.M. Pomirchi, Phys. Sol. State 35, 
1248 (1993). 

23 K. Binder, Z. Phys. B 43,119 (1981). 

24 V. Privman, ed., Finite Size Scaling and Numerical Simu- 
lation of Statistical Systems (World Scientific Publishing, 
Singapore, 1990). 

25 U. Wolff, Nucl. Phys. B 300,501 (1988). 

26 U. Wolff, Phys. Rev. Lett. 62,361 (1989). 

27 C. Kawabata, M. Takeuchi, and AR. Bishop, J. Magn. 
Magn. Mater. 54-57, Pt.2, 871 (1986). 

28 C. Kawabata, M. Takeuchi, and AR. Bishop, J. Magn. 
Magn. Mater. 54-57, Pt.2, 871 (1986). 

29 CM. Wysin and AR. Bishop, Phys. Rev. B 42, 810 
(1990). 

30 M.E. Gouvea and CM. Wysin, Phys. Rev. B 56, 14192 
(1997). 

31 DP. Landau, J. Magn. Magn. Mater. 200,231 (1999). 

32 One MC step includes attempting one single spin move per 
site and one over-relaxation move per site, and individual 
cluster generation until the total number of sites in all 
clusters formed has reached (1/4)L 2 sites. 



